#include "BrayCurtis.h"

BrayCurtis::BrayCurtis()
{
}

BrayCurtis::~BrayCurtis()
{
}

double** BrayCurtis::calc(map<string,double*>& Abundance, int SampleNum, double** BrayCurtisMatrix)
{
    int tPairs=SampleNum*(SampleNum-1)/2;
    double* tX=new double[tPairs];         //sum of abundance of first community
    memset(tX,0,sizeof(double)*tPairs);  
    double* tY=new double[tPairs];         //sum of abundance of second community
    memset(tY,0,sizeof(double)*tPairs);
    double* tC=new double[tPairs];
    memset(tC,0,sizeof(double)*tPairs);
    map<string,double*>::iterator it=Abundance.begin();
    while(it!=Abundance.end())
    {
        for(int i=0;i<SampleNum-1;i++)
            for(int j=i+1;j<SampleNum;j++)
            {
                int tIndex=i*SampleNum-i*(i+1)/2+j-i-1;
                if((*it).second[i]!=0||(*it).second[j]!=0)
                {
                    tX[tIndex]+=(*it).second[i];
                    tY[tIndex]+=(*it).second[j];
                    tC[tIndex]+=((*it).second[i]<(*it).second[j])? (*it).second[i]:(*it).second[j];
                }
            }
        it++;
    }
    for(int i=0;i<SampleNum;i++)
        for(int j=0;j<SampleNum;j++)
        {
            if(i==j)
                BrayCurtisMatrix[i][j]=0;
            if(i>j)
                BrayCurtisMatrix[i][j]=BrayCurtisMatrix[j][i];
            if(i<j)
            {
                int tIndex=i*SampleNum-i*(i+1)/2+j-i-1;
                BrayCurtisMatrix[i][j]=(tX[tIndex]+tY[tIndex]-2.0*tC[tIndex])/(tX[tIndex]+tY[tIndex]);
            }
            //cout<<BrayCurtisMatrix[i][j];
            //if(j==SampleNum-1)
                //cout<<"\n";
            //else
                //cout<<"\t";
        }
    delete[] tX;
    delete[] tY;
    delete[] tC;
    return BrayCurtisMatrix;         
}
